Bee Respiration for Susie and Stacey
Callin Switzer
11/4/2016, Update: 12/13/2016, Update: 28 Jan 2017
Update 28 Jan 2017
Summary of findings:
- % load and total mass are too collinear to include in the same model.
- Treatment and % load are too collinear to include in the same model.
- For the remainder of the modeling, I’ll use only % load
Differences between model with Treatment vs. model with % load
Treatment order is not a significant predictor of wingbeat frequency in the model with percent load (though it is very close, p ~ 0.07)
percent load and treatment suggest very similar models in terms of wingbeat frequency
Treatment order is not a strong predictor of stroke amplitude in the model with percent load
wingbeat freq = size Met rate = mass carrying amplitude = percent loading
Susie’s (hypotheses): lm( frequency ~ ITspan + trt order)
lm( metabolism ~ total mass)
lm( stroke amplitude ~ %load + trt) (those are probably the same thing)
- Redo the model with %load and/or total mass (probably the best predictor of metabolism)
- Check VIF
- Think about how to make figures
How heavy the bee is loaded overall, affects how much the amplitude and frequency change
Questions to look into: ** A really heavy loaded bee will have a smaller change in metabolism vs. light bee ** Suggests a non-linear relationship or random slopes that depend on bee size — the cost for additional weight goes down as you get heavier loaded. — not the case for the change in amplitude (change in load vs. change in amplitude is linear) — not the case for metabolism or wingbeat frequency.
** Hypotheses: 1. Wingbeat freq moves around for a number of reasons (tired, cold, added mass) 2. Stroke amplitude rescues wingbeat frequency 3. As you get really close to your max, adding more weight doesn’t seem to cost as much.
Analyses of respirometry data:
- How does carrying a load during flight affect bumblebees’?
- How does load affect respiratory rate?
- How does load affect wingbeat frequency?
- How does load affect stroke amplitude?
- How does load affect wing velocity?
- Which flight measurements are most closely associated with respiratory rate in bees – wingbeat frequency, stroke amplitude, or wing velocity?
Setup
Install required packages and read in data Define custom function for evaluating VIF with multilevel models
ipak <- function(pkg){
new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
if(length(new.pkg)) install.packages(new.pkg, dependencies = TRUE)
sapply(pkg, require, character.only = TRUE)
}
packages <- c('ggplot2', 'car', 'lme4', 'gsheet', "MASS")
ipak(packages)## ggplot2 car lme4 gsheet MASS
## TRUE TRUE TRUE TRUE TRUEtheme_set(theme_bw())
# read in data -- google sheet called "Bumble mumble grumble"
bdta <- gsheet2tbl("https://docs.google.com/spreadsheets/d/1GUUoFq41Ep3sJNFXiyMcp7fZhYmQCzTcLnhVXdr4WRo/edit?usp=sharing")Function for calculating vif for lmer
(used later)
vif.mer <- function (fit) {
## adapted from rms::vif
v <- vcov(fit)
nam <- names(fixef(fit))
## exclude intercepts
ns <- sum(1 * (nam == "Intercept" | nam == "(Intercept)"))
if (ns > 0) {
v <- v[-(1:ns), -(1:ns), drop = FALSE]
nam <- nam[-(1:ns)]
}
d <- diag(v)^0.5
v <- diag(solve(v/(d %o% d)))
names(v) <- nam
v
}Data Overview
# calculate percent load
bdta$percLoad <- with(bdta, (Mstarved + load) / Mstarved)
# calculate total mass
bdta$totalMass <- with(bdta, Mstarved + load)
summary(bdta)## Bee.ID Treatment.order Treatment Mstarved
## Length:60 Min. :1.0 Length:60 Min. :0.0767
## Class :character 1st Qu.:1.0 Class :character 1st Qu.:0.1213
## Mode :character Median :1.5 Mode :character Median :0.1378
## Mean :1.5 Mean :0.1411
## 3rd Qu.:2.0 3rd Qu.:0.1670
## Max. :2.0 Max. :0.1909
## load IT.Span single.wing.area av.resp..CO2.mL.hr.
## Min. :0.01650 Min. :3.860 Min. :20.03 Min. : 5.507
## 1st Qu.:0.04713 1st Qu.:4.400 1st Qu.:27.93 1st Qu.: 8.935
## Median :0.07840 Median :4.620 Median :30.67 Median :10.666
## Mean :0.08332 Mean :4.594 Mean :30.70 Mean :10.600
## 3rd Qu.:0.11800 3rd Qu.:4.870 3rd Qu.:33.75 3rd Qu.:12.297
## Max. :0.17460 Max. :5.180 Max. :39.05 Max. :16.475
## wbf.aud. stroke.amplitude wing.velocity percLoad
## Min. :163.9 Min. : 96.88 Min. :1.208 Min. :1.143
## 1st Qu.:173.0 1st Qu.:114.16 1st Qu.:1.351 1st Qu.:1.321
## Median :179.1 Median :125.96 Median :1.427 Median :1.576
## Mean :178.9 Mean :123.31 Mean :1.465 Mean :1.600
## 3rd Qu.:185.3 3rd Qu.:133.61 3rd Qu.:1.588 3rd Qu.:1.863
## Max. :194.7 Max. :144.78 Max. :1.843 Max. :2.180
## totalMass
## Min. :0.1062
## 1st Qu.:0.1808
## Median :0.2191
## Mean :0.2245
## 3rd Qu.:0.2646
## Max. :0.3550ggplot(bdta, aes(x = percLoad, fill = Treatment)) +
geom_histogram() ## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# seems like percent load and trt are very related# visualize % load vs total mass
ggplot(bdta, aes(x = percLoad, y = Mstarved + load)) +
geom_point()
# seem to be strongly associated# see number of observations per bee
table(bdta$Bee.ID) # each bee has two observations##
## E01 E03 E05 E09 E10 E11 E12 E16 E20 E24 E26 E28 E30 E31 E32 E33 E35 E36
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## E37 E39 E41 E42 E44 E45 E49 E50 E53 E54 E56 E59
## 2 2 2 2 2 2 2 2 2 2 2 2# get number of unque bees
length(table(bdta$Bee.ID))## [1] 30# get treatment orders
trtOrders_loadedSecond <- sapply(X = unique(bdta$Bee.ID), FUN = function(x){
tmp = bdta[bdta$Bee.ID == x, c("Treatment.order", "Treatment")]
if("2_L" %in% paste(tmp[,1], tmp[,2], sep = "_")){
loadedSecond = TRUE
}
else loadedSecond = FALSE
return(loadedSecond)
})
table(trtOrders_loadedSecond)## trtOrders_loadedSecond
## FALSE
## 30# visualize scatterplot
car::scatterplotMatrix(bdta[, 4:13])
Variance Inflation Factor and PCA
# VIF is high among the bee size predictors
# note: this function is almost the same as car::vif
vif.mer(lmer(av.resp..CO2.mL.hr. ~ Treatment + Treatment.order + Mstarved + IT.Span + totalMass + single.wing.area + percLoad + (1|Bee.ID), data = bdta))## TreatmentUL Treatment.order Mstarved IT.Span
## 16.591645 1.270430 23.112962 6.658432
## totalMass single.wing.area percLoad
## 32.931859 10.903672 38.415618car::scatterplotMatrix(bdta[, c("Mstarved", "IT.Span", "single.wing.area")])
# principle components
aa = prcomp(bdta[, c("Mstarved", "IT.Span", "single.wing.area")], center = TRUE, scale = TRUE)
summary(aa) # 1st pc explains ~95% of the variance in the three measurements of size## Importance of components:
## PC1 PC2 PC3
## Standard deviation 1.6852 0.33843 0.21323
## Proportion of Variance 0.9467 0.03818 0.01516
## Cumulative Proportion 0.9467 0.98484 1.00000biplot(aa) # shows that all three size measurement are correlated
# note, I changed the signs of the predictions so that higher PC1 values
# correspond to bigger bees
p1 = -predict(aa)[,1]
# add PC1 scores to dataset
bdta$size_pca1 = p1
# show scatterplot matrix to see correlations among size predictors
car::scatterplotMatrix(bdta[, c("Mstarved", "IT.Span", "single.wing.area", "size_pca1", "percLoad")])
# check VIF one more time
# VIF is high among the bee size predictors
vif.mer(lmer(av.resp..CO2.mL.hr. ~ Treatment + Treatment.order + size_pca1 + percLoad + totalMass + (1|Bee.ID), data = bdta))## TreatmentUL Treatment.order size_pca1 percLoad
## 16.643643 1.258990 7.448179 30.631628
## totalMass
## 26.790114# remove treatment
vif.mer(lmer(av.resp..CO2.mL.hr. ~ Treatment.order + size_pca1 + percLoad + totalMass + (1|Bee.ID), data = bdta))## Treatment.order size_pca1 percLoad totalMass
## 1.014834 7.698050 20.910198 25.813622# looks like we cannot have total mass and percentLoad in the same model, according to VIF
# remove total mass
vif.mer(lmer(av.resp..CO2.mL.hr. ~ Treatment.order + size_pca1 + percLoad + (1|Bee.ID), data = bdta))## Treatment.order size_pca1 percLoad
## 1.003811 1.004814 1.008624Multilevel models
# a few house-keeping issues to make models easier to read
# convert trt order to factor
bdta$Treatment.order <- as.factor(as.character(bdta$Treatment.order))
# change reference level to unloaded
bdta$Treatment <- factor(bdta$Treatment, levels = c("UL", "L"))Resp Rate
# make a full model with all two-way interactions
m1 <- lmer(av.resp..CO2.mL.hr. ~ (size_pca1 + Treatment.order + percLoad)^2 + (1|Bee.ID), data = bdta)
summary(m1)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## av.resp..CO2.mL.hr. ~ (size_pca1 + Treatment.order + percLoad)^2 +
## (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 196.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.5230 -0.4541 -0.0409 0.4813 2.2994
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.2498 1.1179
## Residual 0.7711 0.8781
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.0096 1.3291 2.264
## size_pca1 0.6926 0.4455 1.555
## Treatment.order2 1.9553 2.0573 0.950
## percLoad 4.8854 0.8030 6.084
## size_pca1:Treatment.order2 -0.0830 0.1451 -0.572
## size_pca1:percLoad 0.3533 0.2491 1.418
## Treatment.order2:percLoad -1.4903 1.2730 -1.171
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 percLd s_1:T. sz_1:L
## size_pca1 -0.223
## Trtmnt.rdr2 -0.860 0.246
## percLoad -0.980 0.204 0.866
## sz_pc1:Tr.2 0.297 -0.379 -0.273 -0.297
## sz_pc1:prcL 0.148 -0.933 -0.184 -0.125 0.228
## Trtmnt.r2:L 0.853 -0.233 -0.994 -0.868 0.271 0.169### LRT's for interactions
m2.0 <- update(m1, .~. - Treatment.order:percLoad)
anova(m1, m2.0) # likelihood ratio test for interaction of treatment order## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m2.0: av.resp..CO2.mL.hr. ~ size_pca1 + Treatment.order + percLoad +
## m2.0: (1 | Bee.ID) + size_pca1:Treatment.order + size_pca1:percLoad
## m1: av.resp..CO2.mL.hr. ~ (size_pca1 + Treatment.order + percLoad)^2 +
## m1: (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.0 8 208.02 224.78 -96.011 192.02
## m1 9 208.55 227.40 -95.275 190.55 1.472 1 0.225summary(m2.0)## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ size_pca1 + Treatment.order + percLoad +
## (1 | Bee.ID) + size_pca1:Treatment.order + size_pca1:percLoad
## Data: bdta
##
## REML criterion at convergence: 199.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.61501 -0.46662 -0.04214 0.53248 2.31063
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.3003 1.1403
## Residual 0.7595 0.8715
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 4.33950 0.69115 6.279
## size_pca1 0.57136 0.43113 1.325
## Treatment.order2 -0.43829 0.22783 -1.924
## percLoad 4.06743 0.39578 10.277
## size_pca1:Treatment.order2 -0.03698 0.13862 -0.267
## size_pca1:percLoad 0.40235 0.24379 1.650
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 percLd s_1:T.
## size_pca1 -0.048
## Trtmnt.rdr2 -0.219 0.137
## percLoad -0.924 0.003 0.055
## sz_pc1:Tr.2 0.131 -0.337 -0.036 -0.131
## sz_pc1:prcL 0.007 -0.930 -0.144 0.045 0.192m2.1 <- update(m2.0,.~. - size_pca1:Treatment.order)
anova(m2.0, m2.1)## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m2.1: av.resp..CO2.mL.hr. ~ size_pca1 + Treatment.order + percLoad +
## m2.1: (1 | Bee.ID) + size_pca1:percLoad
## m2.0: av.resp..CO2.mL.hr. ~ size_pca1 + Treatment.order + percLoad +
## m2.0: (1 | Bee.ID) + size_pca1:Treatment.order + size_pca1:percLoad
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.1 7 206.10 220.76 -96.051 192.10
## m2.0 8 208.02 224.78 -96.011 192.02 0.0814 1 0.7754summary(m2.1)## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ size_pca1 + Treatment.order + percLoad +
## (1 | Bee.ID) + size_pca1:percLoad
## Data: bdta
##
## REML criterion at convergence: 197.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.68107 -0.47084 -0.05908 0.51796 2.31287
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.3135 1.1461
## Residual 0.7334 0.8564
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 4.3659 0.6752 6.466
## size_pca1 0.5329 0.4000 1.332
## Treatment.order2 -0.4405 0.2237 -1.969
## percLoad 4.0521 0.3858 10.502
## size_pca1:percLoad 0.4146 0.2352 1.763
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 percLd
## size_pca1 -0.004
## Trtmnt.rdr2 -0.216 0.132
## percLoad -0.921 -0.044 0.051
## sz_pc1:prcL -0.019 -0.935 -0.140 0.072m2.2 <- update(m2.1, .~. - size_pca1:percLoad)
anova(m2.1, m2.2) ## drop all interactions## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m2.2: av.resp..CO2.mL.hr. ~ size_pca1 + Treatment.order + percLoad +
## m2.2: (1 | Bee.ID)
## m2.1: av.resp..CO2.mL.hr. ~ size_pca1 + Treatment.order + percLoad +
## m2.1: (1 | Bee.ID) + size_pca1:percLoad
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.2 6 207.38 219.94 -97.688 195.38
## m2.1 7 206.10 220.76 -96.051 192.10 3.2723 1 0.07046 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m2.2)## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ size_pca1 + Treatment.order + percLoad +
## (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 199.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.44974 -0.62601 -0.00733 0.59569 2.08255
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.296 1.138
## Residual 0.785 0.886
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 4.3843 0.6950 6.309
## size_pca1 1.1922 0.1423 8.378
## Treatment.order2 -0.3854 0.2292 -1.681
## percLoad 4.0055 0.3976 10.074
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2
## size_pca1 -0.064
## Trtmnt.rdr2 -0.221 0.004
## percLoad -0.925 0.069 0.062# renaming model to simplify later typing
m2 <- m2.2
##### LRTs for main effects
## Treatment Order
m3 <- update(m2, .~. - Treatment.order)
anova(m2, m3, test = "Chi") # drop trt order (different than model with treatment)## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m3: av.resp..CO2.mL.hr. ~ size_pca1 + percLoad + (1 | Bee.ID)
## m2: av.resp..CO2.mL.hr. ~ size_pca1 + Treatment.order + percLoad +
## m2: (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m3 5 208.25 218.73 -99.127 198.25
## m2 6 207.38 219.94 -97.688 195.38 2.8789 1 0.08975 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# LRT for size
m4 <- update(m3, .~. - size_pca1)
anova(m3, m4) # keep size## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m4: av.resp..CO2.mL.hr. ~ percLoad + (1 | Bee.ID)
## m3: av.resp..CO2.mL.hr. ~ size_pca1 + percLoad + (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m4 4 244.02 252.40 -118.011 236.02
## m3 5 208.25 218.73 -99.127 198.25 37.767 1 7.97e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# LRT for Treatment (load)
m5 <- update(m3, .~. - percLoad)
anova(m3, m5) # keep percent load## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m5: av.resp..CO2.mL.hr. ~ size_pca1 + (1 | Bee.ID)
## m3: av.resp..CO2.mL.hr. ~ size_pca1 + percLoad + (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m5 4 255.27 263.64 -123.634 247.27
## m3 5 208.25 218.73 -99.127 198.25 49.013 1 2.543e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# summarize final model for paper
summary(m3)## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ size_pca1 + percLoad + (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 201.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.33113 -0.59393 0.02882 0.58497 2.24345
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.2695 1.1267
## Residual 0.8349 0.9137
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 4.1216 0.6955 5.926
## size_pca1 1.1932 0.1423 8.388
## percLoad 4.0493 0.4087 9.907
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1
## size_pca1 -0.067
## percLoad -0.940 0.071# write output
summary(m3)$coefficients ## Estimate Std. Error t value
## (Intercept) 4.121569 0.6955406 5.925706
## size_pca1 1.193243 0.1422611 8.387692
## percLoad 4.049269 0.4087169 9.907270write.csv( summary(m3)$coefficients, file = "RespCoefs_percLoad.csv" )
# diagnostics
# qq plot
qqnorm(resid(m3), main = "")
qqline(resid(m3)) # good
# residual plot
plot(fitted(m3), resid(m3), xlab = "fitted", ylab = "residuals")
abline(0,0)
# QQPlot for group-level effects
qqnorm(ranef(m3)$Bee.ID[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(m3)$Bee.ID[[1]]) # looks good
Wingbeat Freq
# make a full model with all two-way interactions
m1 <- lmer(wbf.aud. ~ (size_pca1 + Treatment.order + percLoad)^2 + (1|Bee.ID), data = bdta)
summary(m1)## Linear mixed model fit by REML ['lmerMod']
## Formula: wbf.aud. ~ (size_pca1 + Treatment.order + percLoad)^2 + (1 |
## Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 349.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.80184 -0.60216 -0.06242 0.61714 1.64271
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 10.18 3.190
## Residual 20.62 4.541
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 173.7585 5.7634 30.149
## size_pca1 -4.8135 2.1942 -2.194
## Treatment.order2 -5.8644 8.3595 -0.702
## percLoad 4.3885 3.4953 1.256
## size_pca1:Treatment.order2 0.9035 0.7397 1.221
## size_pca1:percLoad 0.9858 1.2643 0.780
## Treatment.order2:percLoad 1.4493 5.1496 0.281
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 percLd s_1:T. sz_1:L
## size_pca1 -0.198
## Trtmnt.rdr2 -0.818 0.226
## percLoad -0.984 0.172 0.814
## sz_pc1:Tr.2 0.254 -0.376 -0.223 -0.252
## sz_pc1:prcL 0.127 -0.958 -0.171 -0.098 0.216
## Trtmnt.r2:L 0.807 -0.209 -0.990 -0.818 0.221 0.152### LRT's for interactions
m2.0 <- update(m1, .~. - size_pca1:percLoad )
anova(m1, m2.0)## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m2.0: wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID) +
## m2.0: size_pca1:Treatment.order + Treatment.order:percLoad
## m1: wbf.aud. ~ (size_pca1 + Treatment.order + percLoad)^2 + (1 |
## m1: Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.0 8 381.69 398.44 -182.84 365.69
## m1 9 382.99 401.83 -182.49 364.99 0.7005 1 0.4026summary(m2.0)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID) +
## size_pca1:Treatment.order + Treatment.order:percLoad
## Data: bdta
##
## REML criterion at convergence: 352.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.72938 -0.54168 0.01881 0.54877 1.67034
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 10.08 3.175
## Residual 20.49 4.526
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 173.1944 5.6964 30.404
## size_pca1 -3.1748 0.6277 -5.058
## Treatment.order2 -4.7581 8.2056 -0.580
## percLoad 4.6533 3.4658 1.343
## size_pca1:Treatment.order2 0.7793 0.7199 1.082
## Treatment.order2:percLoad 0.8413 5.0701 0.166
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 percLd s_1:T.
## size_pca1 -0.268
## Trtmnt.rdr2 -0.815 0.221
## percLoad -0.984 0.272 0.813
## sz_pc1:Tr.2 0.234 -0.606 -0.194 -0.238
## Trtmnt.r2:L 0.803 -0.222 -0.990 -0.816 0.194m2.1 <- update(m2.0,.~. - size_pca1:Treatment.order)
anova(m2.0, m2.1)## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m2.1: wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID) +
## m2.1: Treatment.order:percLoad
## m2.0: wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID) +
## m2.0: size_pca1:Treatment.order + Treatment.order:percLoad
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.1 7 380.97 395.64 -183.49 366.97
## m2.0 8 381.69 398.44 -182.84 365.69 1.2893 1 0.2562summary(m2.1)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID) +
## Treatment.order:percLoad
## Data: bdta
##
## REML criterion at convergence: 354.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.60845 -0.62485 0.06463 0.64370 1.69662
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 10.35 3.218
## Residual 20.39 4.515
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 171.6958 5.5454 30.962
## size_pca1 -2.7623 0.5021 -5.501
## Treatment.order2 -2.9560 8.0745 -0.366
## percLoad 5.5797 3.3706 1.655
## Treatment.order2:percLoad -0.2770 4.9889 -0.056
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 percLd
## size_pca1 -0.163
## Trtmnt.rdr2 -0.808 0.133
## percLoad -0.983 0.165 0.806
## Trtmnt.r2:L 0.796 -0.133 -0.989 -0.809m2.2 <- update(m2.1, .~. - Treatment.order:percLoad)
anova(m2.1, m2.2) ## drop all interactions## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m2.2: wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID)
## m2.1: wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID) +
## m2.1: Treatment.order:percLoad
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.2 6 378.98 391.54 -183.49 366.98
## m2.1 7 380.97 395.64 -183.49 366.97 0.0018 1 0.9666summary(m2.2)## Linear mixed model fit by REML ['lmerMod']
## Formula: wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 360
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.64247 -0.64283 0.06575 0.64332 1.72090
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 9.874 3.142
## Residual 20.221 4.497
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 171.9557 3.3390 51.50
## size_pca1 -2.7663 0.4908 -5.64
## Treatment.order2 -3.3999 1.1632 -2.92
## percLoad 5.4191 1.9692 2.75
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2
## size_pca1 -0.095
## Trtmnt.rdr2 -0.231 0.006
## percLoad -0.954 0.099 0.060# renaming model to simplify later typing
m2 <- m2.2
##### LRTs for main effects
## Treatment Order
m3 <- update(m2, .~. - Treatment.order)
anova(m2, m3, test = "Chi")## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m3: wbf.aud. ~ size_pca1 + percLoad + (1 | Bee.ID)
## m2: wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m3 5 384.96 395.43 -187.48 374.96
## m2 6 378.98 391.54 -183.49 366.98 7.9842 1 0.004719 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# LRT for size
m4 <- update(m2, .~. - size_pca1)
anova(m2, m4)## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m4: wbf.aud. ~ Treatment.order + percLoad + (1 | Bee.ID)
## m2: wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m4 5 399.76 410.23 -194.88 389.76
## m2 6 378.98 391.54 -183.49 366.98 22.783 1 1.813e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# LRT for Treatment (load)
m5 <- update(m2, .~. - percLoad)
anova(m2, m5)## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m5: wbf.aud. ~ size_pca1 + Treatment.order + (1 | Bee.ID)
## m2: wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m5 5 383.97 394.45 -186.99 373.97
## m2 6 378.98 391.54 -183.49 366.98 6.9967 1 0.008166 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# summarize final model for paper
summary(m2)## Linear mixed model fit by REML ['lmerMod']
## Formula: wbf.aud. ~ size_pca1 + Treatment.order + percLoad + (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 360
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.64247 -0.64283 0.06575 0.64332 1.72090
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 9.874 3.142
## Residual 20.221 4.497
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 171.9557 3.3390 51.50
## size_pca1 -2.7663 0.4908 -5.64
## Treatment.order2 -3.3999 1.1632 -2.92
## percLoad 5.4191 1.9692 2.75
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2
## size_pca1 -0.095
## Trtmnt.rdr2 -0.231 0.006
## percLoad -0.954 0.099 0.060# write output
summary(m2)$coefficients ## Estimate Std. Error t value
## (Intercept) 171.955651 3.3390335 51.498630
## size_pca1 -2.766273 0.4908262 -5.635952
## Treatment.order2 -3.399880 1.1631609 -2.922966
## percLoad 5.419093 1.9691779 2.751957write.csv( summary(m2)$coefficients, file = "FreqCoefs_percLoad.csv" )
# diagnostics
# qq plot
qqnorm(resid(m2), main = "")
qqline(resid(m2)) # good
# residual plot
plot(fitted(m2), resid(m2), xlab = "fitted", ylab = "residuals")
abline(0,0)
# QQPlot for group-level effects
qqnorm(ranef(m2)$Bee.ID[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(m2)$Bee.ID[[1]]) # looks good
Stroke Amplitude
# make a full model with all two-way interactions
m1 <- lmer(stroke.amplitude ~ (size_pca1 + Treatment.order + percLoad)^2 + (1|Bee.ID), data = bdta)
summary(m1)## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ (size_pca1 + Treatment.order + percLoad)^2 +
## (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 359.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.89733 -0.52961 -0.01694 0.40969 1.99945
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 14.12 3.758
## Residual 23.72 4.870
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 61.9443 6.3305 9.785
## size_pca1 6.0503 2.3672 2.556
## Treatment.order2 7.4111 9.2790 0.799
## percLoad 38.4339 3.8380 10.014
## size_pca1:Treatment.order2 -1.1440 0.7947 -1.440
## size_pca1:percLoad -2.6118 1.3601 -1.920
## Treatment.order2:percLoad -5.0714 5.7204 -0.887
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 percLd s_1:T. sz_1:L
## size_pca1 -0.201
## Trtmnt.rdr2 -0.825 0.229
## percLoad -0.984 0.176 0.822
## sz_pc1:Tr.2 0.259 -0.377 -0.230 -0.258
## sz_pc1:prcL 0.129 -0.956 -0.172 -0.102 0.217
## Trtmnt.r2:L 0.814 -0.212 -0.991 -0.826 0.227 0.154### LRT's for interactions
m2.0 <- update(m1, .~. - Treatment.order:percLoad )
anova(m1, m2.0)## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m2.0: stroke.amplitude ~ size_pca1 + Treatment.order + percLoad + (1 |
## m2.0: Bee.ID) + size_pca1:Treatment.order + size_pca1:percLoad
## m1: stroke.amplitude ~ (size_pca1 + Treatment.order + percLoad)^2 +
## m1: (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.0 8 393.14 409.89 -188.57 377.14
## m1 9 394.31 413.16 -188.16 376.31 0.8251 1 0.3637summary(m2.0)## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ size_pca1 + Treatment.order + percLoad + (1 |
## Bee.ID) + size_pca1:Treatment.order + size_pca1:percLoad
## Data: bdta
##
## REML criterion at convergence: 365.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.89468 -0.52361 -0.02085 0.43365 2.07383
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 14.75 3.841
## Residual 23.18 4.814
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 66.4391 3.6420 18.242
## size_pca1 5.5979 2.2913 2.443
## Treatment.order2 -0.7365 1.2581 -0.585
## percLoad 35.6715 2.1438 16.639
## size_pca1:Treatment.order2 -0.9857 0.7651 -1.288
## size_pca1:percLoad -2.4201 1.3296 -1.820
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 percLd s_1:T.
## size_pca1 -0.050
## Trtmnt.rdr2 -0.227 0.138
## percLoad -0.950 0.002 0.054
## sz_pc1:Tr.2 0.132 -0.345 -0.035 -0.128
## sz_pc1:prcL 0.006 -0.955 -0.142 0.046 0.189m2.1 <- update(m2.0,.~. - size_pca1:Treatment.order)
anova(m2.0, m2.1)## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m2.1: stroke.amplitude ~ size_pca1 + Treatment.order + percLoad + (1 |
## m2.1: Bee.ID) + size_pca1:percLoad
## m2.0: stroke.amplitude ~ size_pca1 + Treatment.order + percLoad + (1 |
## m2.0: Bee.ID) + size_pca1:Treatment.order + size_pca1:percLoad
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.1 7 392.99 407.65 -189.49 378.99
## m2.0 8 393.14 409.89 -188.57 377.14 1.8488 1 0.1739summary(m2.1)## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ size_pca1 + Treatment.order + percLoad + (1 |
## Bee.ID) + size_pca1:percLoad
## Data: bdta
##
## REML criterion at convergence: 368.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.92050 -0.50058 -0.06054 0.48999 2.08830
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 14.24 3.773
## Residual 23.88 4.887
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 67.1284 3.6558 18.362
## size_pca1 4.5889 2.1787 2.106
## Treatment.order2 -0.7942 1.2763 -0.622
## percLoad 35.2724 2.1552 16.367
## size_pca1:percLoad -2.1023 1.3241 -1.588
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 percLd
## size_pca1 -0.005
## Trtmnt.rdr2 -0.225 0.135
## percLoad -0.950 -0.046 0.050
## sz_pc1:prcL -0.019 -0.966 -0.138 0.072m2.2 <- update(m2.1, .~. - size_pca1:percLoad)
anova(m2.1, m2.2) ## drop all interactions## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m2.2: stroke.amplitude ~ size_pca1 + Treatment.order + percLoad + (1 |
## m2.2: Bee.ID)
## m2.1: stroke.amplitude ~ size_pca1 + Treatment.order + percLoad + (1 |
## m2.1: Bee.ID) + size_pca1:percLoad
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.2 6 393.68 406.24 -190.84 381.68
## m2.1 7 392.99 407.65 -189.49 378.99 2.6909 1 0.1009summary(m2.2)## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ size_pca1 + Treatment.order + percLoad + (1 |
## Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 373.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.87059 -0.50457 -0.03627 0.58077 2.02395
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 14.31 3.782
## Residual 24.73 4.973
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 67.0491 3.7152 18.047
## size_pca1 1.2463 0.5668 2.199
## Treatment.order2 -1.0742 1.2864 -0.835
## percLoad 35.5004 2.1858 16.241
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2
## size_pca1 -0.091
## Trtmnt.rdr2 -0.230 0.006
## percLoad -0.952 0.096 0.060# renaming model to simplify later typing
m2 <- m2.2
##### LRTs for main effects
## size
m3 <- update(m2, .~. - size_pca1)
anova(m2, m3, test = "Chi") # keep size_pca1## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m3: stroke.amplitude ~ Treatment.order + percLoad + (1 | Bee.ID)
## m2: stroke.amplitude ~ size_pca1 + Treatment.order + percLoad + (1 |
## m2: Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m3 5 396.50 406.97 -193.25 386.50
## m2 6 393.68 406.24 -190.84 381.68 4.8203 1 0.02813 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m3)## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ Treatment.order + percLoad + (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 379
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.79024 -0.53442 0.00173 0.52524 1.85491
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 17.63 4.199
## Residual 24.87 4.987
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 67.501 3.741 18.04
## Treatment.order2 -1.084 1.290 -0.84
## percLoad 35.221 2.193 16.06
##
## Correlation of Fixed Effects:
## (Intr) Trtm.2
## Trtmnt.rdr2 -0.229
## percLoad -0.948 0.060# trt order
m4 <- update(m2, .~. - Treatment.order)
anova(m2, m4) # drop treatment order## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m4: stroke.amplitude ~ size_pca1 + percLoad + (1 | Bee.ID)
## m2: stroke.amplitude ~ size_pca1 + Treatment.order + percLoad + (1 |
## m2: Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m4 5 392.42 402.89 -191.21 382.42
## m2 6 393.68 406.24 -190.84 381.68 0.7374 1 0.3905# LRT for Treatment (load)
m5 <- update(m4, .~. - percLoad)
anova(m4, m5)## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m5: stroke.amplitude ~ size_pca1 + (1 | Bee.ID)
## m4: stroke.amplitude ~ size_pca1 + percLoad + (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m5 4 478.07 486.45 -235.04 470.07
## m4 5 392.42 402.89 -191.21 382.42 87.657 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# summarize final model for paper
summary(m4)## Linear mixed model fit by REML ['lmerMod']
## Formula: stroke.amplitude ~ size_pca1 + percLoad + (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 376.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.96941 -0.48117 -0.05909 0.48300 1.94459
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 14.53 3.811
## Residual 24.42 4.941
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 66.3057 3.5966 18.436
## size_pca1 1.2495 0.5675 2.202
## percLoad 35.6293 2.1693 16.424
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1
## size_pca1 -0.091
## percLoad -0.965 0.095# write output
summary(m4)$coefficients ## Estimate Std. Error t value
## (Intercept) 66.305690 3.5965822 18.435750
## size_pca1 1.249463 0.5674547 2.201872
## percLoad 35.629329 2.1693080 16.424283write.csv( summary(m4)$coefficients, file = "AmpCoefs_percLoad.csv" )
# diagnostics
# qq plot
qqnorm(resid(m4), main = "")
qqline(resid(m4)) # ok
# residual plot
plot(fitted(m4), resid(m3), xlab = "fitted", ylab = "residuals")
abline(0,0)
# QQPlot for group-level effects
qqnorm(ranef(m4)$Bee.ID[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(m4)$Bee.ID[[1]]) # looks good
Wing Velocity
# make a full model with all two-way interactions
m1 <- lmer(wing.velocity ~ (size_pca1 + Treatment.order + percLoad)^2 + (1|Bee.ID), data = bdta)
summary(m1)## Linear mixed model fit by REML ['lmerMod']
## Formula: wing.velocity ~ (size_pca1 + Treatment.order + percLoad)^2 +
## (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: -78.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.00651 -0.41257 -0.06598 0.54758 2.20129
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 0.003893 0.06239
## Residual 0.005935 0.07704
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.14983 0.10145 21.190
## size_pca1 -0.12816 0.03757 -3.411
## Treatment.order2 -0.11249 0.14953 -0.752
## percLoad -0.41933 0.06150 -6.819
## size_pca1:Treatment.order2 0.01603 0.01258 1.274
## size_pca1:percLoad 0.05043 0.02155 2.340
## Treatment.order2:percLoad 0.05830 0.09222 0.632
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 percLd s_1:T. sz_1:L
## size_pca1 -0.203
## Trtmnt.rdr2 -0.828 0.230
## percLoad -0.984 0.178 0.826
## sz_pc1:Tr.2 0.262 -0.377 -0.233 -0.261
## sz_pc1:prcL 0.131 -0.954 -0.173 -0.103 0.218
## Trtmnt.r2:L 0.818 -0.214 -0.991 -0.830 0.231 0.155### LRT's for interactions
m2.0 <- update(m1, .~. - Treatment.order:percLoad )
anova(m1, m2.0)## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m2.0: wing.velocity ~ size_pca1 + Treatment.order + percLoad + (1 |
## m2.0: Bee.ID) + size_pca1:Treatment.order + size_pca1:percLoad
## m1: wing.velocity ~ (size_pca1 + Treatment.order + percLoad)^2 +
## m1: (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.0 8 -103.25 -86.492 59.624 -119.25
## m1 9 -101.66 -82.811 59.830 -119.66 0.4129 1 0.5205summary(m2.0)## Linear mixed model fit by REML ['lmerMod']
## Formula: wing.velocity ~ size_pca1 + Treatment.order + percLoad + (1 |
## Bee.ID) + size_pca1:Treatment.order + size_pca1:percLoad
## Data: bdta
##
## REML criterion at convergence: -80.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.07164 -0.43946 -0.00926 0.53663 2.21665
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 0.003959 0.06292
## Residual 0.005800 0.07616
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.09780 0.05777 36.31
## size_pca1 -0.12306 0.03633 -3.39
## Treatment.order2 -0.01883 0.01990 -0.95
## percLoad -0.38734 0.03397 -11.40
## size_pca1:Treatment.order2 0.01421 0.01210 1.17
## size_pca1:percLoad 0.04829 0.02106 2.29
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 percLd s_1:T.
## size_pca1 -0.050
## Trtmnt.rdr2 -0.227 0.138
## percLoad -0.949 0.002 0.054
## sz_pc1:Tr.2 0.132 -0.345 -0.035 -0.129
## sz_pc1:prcL 0.006 -0.954 -0.142 0.046 0.190m2.1 <- update(m2.0,.~. - size_pca1:Treatment.order)
anova(m2.0, m2.1)## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m2.1: wing.velocity ~ size_pca1 + Treatment.order + percLoad + (1 |
## m2.1: Bee.ID) + size_pca1:percLoad
## m2.0: wing.velocity ~ size_pca1 + Treatment.order + percLoad + (1 |
## m2.0: Bee.ID) + size_pca1:Treatment.order + size_pca1:percLoad
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.1 7 -103.70 -89.044 58.852 -117.70
## m2.0 8 -103.25 -86.492 59.624 -119.25 1.5428 1 0.2142summary(m2.1)## Linear mixed model fit by REML ['lmerMod']
## Formula: wing.velocity ~ size_pca1 + Treatment.order + percLoad + (1 |
## Bee.ID) + size_pca1:percLoad
## Data: bdta
##
## REML criterion at convergence: -86.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.08227 -0.51145 -0.04384 0.65437 2.23563
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 0.003885 0.06233
## Residual 0.005902 0.07682
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.08849 0.05768 36.21
## size_pca1 -0.10840 0.03436 -3.16
## Treatment.order2 -0.01800 0.02006 -0.90
## percLoad -0.38198 0.03395 -11.25
## size_pca1:percLoad 0.04363 0.02084 2.09
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 Trtm.2 percLd
## size_pca1 -0.005
## Trtmnt.rdr2 -0.224 0.135
## percLoad -0.948 -0.046 0.050
## sz_pc1:prcL -0.019 -0.965 -0.138 0.072m2.2 <- update(m2.1, .~. - size_pca1:percLoad)
anova(m2.1, m2.2) ## don't drop size:percLoad interaction## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m2.2: wing.velocity ~ size_pca1 + Treatment.order + percLoad + (1 |
## m2.2: Bee.ID)
## m2.1: wing.velocity ~ size_pca1 + Treatment.order + percLoad + (1 |
## m2.1: Bee.ID) + size_pca1:percLoad
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2.2 6 -101.16 -88.596 56.581 -113.16
## m2.1 7 -103.70 -89.044 58.852 -117.70 4.5421 1 0.03307 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# renaming model to simplify later typing
m2 <- m2.1
##### LRTs for main effects
# trt order
m3 <- update(m2, .~. - Treatment.order)
anova(m2, m3, test = "Chi") # drop trt order## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## m3: wing.velocity ~ size_pca1 + percLoad + (1 | Bee.ID) + size_pca1:percLoad
## m2: wing.velocity ~ size_pca1 + Treatment.order + percLoad + (1 |
## m2: Bee.ID) + size_pca1:percLoad
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m3 6 -104.83 -92.261 58.413 -116.83
## m2 7 -103.70 -89.044 58.852 -117.70 0.8775 1 0.3489summary(m3)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## wing.velocity ~ size_pca1 + percLoad + (1 | Bee.ID) + size_pca1:percLoad
## Data: bdta
##
## REML criterion at convergence: -91.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.18903 -0.53485 -0.00798 0.63629 2.12781
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 0.003898 0.06243
## Residual 0.005865 0.07658
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.07698 0.05606 37.05
## size_pca1 -0.10424 0.03395 -3.07
## percLoad -0.38053 0.03381 -11.25
## size_pca1:percLoad 0.04104 0.02058 1.99
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 percLd
## size_pca1 0.026
## percLoad -0.963 -0.053
## sz_pc1:prcL -0.052 -0.964 0.080# summarize final model for paper
summary(m3)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## wing.velocity ~ size_pca1 + percLoad + (1 | Bee.ID) + size_pca1:percLoad
## Data: bdta
##
## REML criterion at convergence: -91.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.18903 -0.53485 -0.00798 0.63629 2.12781
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 0.003898 0.06243
## Residual 0.005865 0.07658
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.07698 0.05606 37.05
## size_pca1 -0.10424 0.03395 -3.07
## percLoad -0.38053 0.03381 -11.25
## size_pca1:percLoad 0.04104 0.02058 1.99
##
## Correlation of Fixed Effects:
## (Intr) sz_pc1 percLd
## size_pca1 0.026
## percLoad -0.963 -0.053
## sz_pc1:prcL -0.052 -0.964 0.080# write output
summary(m3)$coefficients ## Estimate Std. Error t value
## (Intercept) 2.07697845 0.05606361 37.046819
## size_pca1 -0.10424096 0.03394778 -3.070627
## percLoad -0.38052624 0.03381059 -11.254646
## size_pca1:percLoad 0.04104281 0.02058223 1.994090write.csv( summary(m3)$coefficients, file = "VelocCoefs_percLoad.csv" )
# diagnostics
# qq plot
qqnorm(resid(m3), main = "")
qqline(resid(m3)) # good
# residual plot
plot(fitted(m3), resid(m3), xlab = "fitted", ylab = "residuals")
abline(0,0)
# QQPlot for group-level effects
qqnorm(ranef(m3)$Bee.ID[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(m3)$Bee.ID[[1]]) # looks good
Part II: Which flight measurements are most closely associated with
respiratory rate in bees – wingbeat frequency, stroke amplitude, or
wing velocity?
Multilevel model approach to see the relatedness of the response variables
- Wingbeat frequency
- Stroke amplitude
- Wing velocity
- Resp. Rate
# predicting avg resp rate, including stroke amplitude, wing velocity and wingbeat frequency
# looks like we have to remove wing velocity
vif.mer(lmer(av.resp..CO2.mL.hr. ~ stroke.amplitude + wing.velocity + wbf.aud. + percLoad + Treatment.order + size_pca1 + (1|Bee.ID), data = bdta))## stroke.amplitude wing.velocity wbf.aud. percLoad
## 190.642034 133.592345 11.644262 10.853222
## Treatment.order2 size_pca1
## 1.262757 1.364694vif.mer(lmer(av.resp..CO2.mL.hr. ~ stroke.amplitude + wbf.aud. + percLoad + Treatment.order + size_pca1 + (1|Bee.ID), data = bdta)) # still worryingly high for percLoad and stroke amplitude## stroke.amplitude wbf.aud. percLoad Treatment.order2
## 7.741481 1.830209 8.182520 1.255652
## size_pca1
## 1.355159mm1 <- lmer(av.resp..CO2.mL.hr. ~ (stroke.amplitude + wbf.aud. + percLoad + Treatment.order + size_pca1)^2 + (1|Bee.ID), data = bdta)## Warning: Some predictor variables are on very different scales: consider
## rescalingsummary(mm1) # warning says to rescale## Linear mixed model fit by REML ['lmerMod']
## Formula: av.resp..CO2.mL.hr. ~ (stroke.amplitude + wbf.aud. + percLoad +
## Treatment.order + size_pca1)^2 + (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 211.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.4812 -0.3747 -0.1374 0.2959 2.1473
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 0.9021 0.9498
## Residual 0.4039 0.6356
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -59.720814 44.421604 -1.344
## stroke.amplitude 0.359487 0.608891 0.590
## wbf.aud. 0.278405 0.254698 1.093
## percLoad 11.733723 29.320882 0.400
## Treatment.order2 9.248884 8.744342 1.058
## size_pca1 3.953128 2.787254 1.418
## stroke.amplitude:wbf.aud. -0.001567 0.003489 -0.449
## stroke.amplitude:percLoad -0.089640 0.050249 -1.784
## stroke.amplitude:Treatment.order2 0.019955 0.036748 0.543
## stroke.amplitude:size_pca1 0.010180 0.016490 0.617
## wbf.aud.:percLoad 0.033891 0.149234 0.227
## wbf.aud.:Treatment.order2 -0.047825 0.043855 -1.091
## wbf.aud.:size_pca1 -0.017492 0.013469 -1.299
## percLoad:Treatment.order2 -1.939776 1.533297 -1.265
## percLoad:size_pca1 -0.158919 0.670849 -0.237
## Treatment.order2:size_pca1 -0.455025 0.180000 -2.528##
## Correlation matrix not shown by default, as p = 16 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it## fit warnings:
## Some predictor variables are on very different scales: consider rescaling# rescale numeric predictors
bdta_scaled = scale(bdta[, c("av.resp..CO2.mL.hr.",
"stroke.amplitude",
"wbf.aud.",
"percLoad",
"size_pca1")],
center = TRUE,
scale = TRUE)
colnames(bdta_scaled) <- paste0(colnames(bdta_scaled), "_scaled")
bdta = cbind(bdta, bdta_scaled)
mm1 <- lmer(av.resp..CO2.mL.hr. ~ (stroke.amplitude_scaled + wbf.aud._scaled + percLoad + Treatment.order + size_pca1_scaled)^2 + (1|Bee.ID), data = bdta)
summary(mm1)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## av.resp..CO2.mL.hr. ~ (stroke.amplitude_scaled + wbf.aud._scaled +
## percLoad + Treatment.order + size_pca1_scaled)^2 + (1 | Bee.ID)
## Data: bdta
##
## REML criterion at convergence: 160.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.4812 -0.3747 -0.1374 0.2959 2.1473
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 0.9021 0.9498
## Residual 0.4039 0.6356
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -0.1472 2.1940 -0.067
## stroke.amplitude_scaled 0.9720 0.9560 1.017
## wbf.aud._scaled 0.6494 1.7496 0.371
## percLoad 6.7446 1.4013 4.813
## Treatment.order2 3.1524 2.4584 1.282
## size_pca1_scaled 3.5029 1.7806 1.967
## stroke.amplitude_scaled:wbf.aud._scaled -0.1466 0.3266 -0.449
## stroke.amplitude_scaled:percLoad -1.1008 0.6171 -1.784
## stroke.amplitude_scaled:Treatment.order2 0.2451 0.4513 0.543
## stroke.amplitude_scaled:size_pca1_scaled 0.2107 0.3413 0.617
## wbf.aud._scaled:percLoad 0.2583 1.1374 0.227
## wbf.aud._scaled:Treatment.order2 -0.3645 0.3342 -1.091
## wbf.aud._scaled:size_pca1_scaled -0.2247 0.1730 -1.299
## percLoad:Treatment.order2 -1.9398 1.5333 -1.265
## percLoad:size_pca1_scaled -0.2678 1.1305 -0.237
## Treatment.order2:size_pca1_scaled -0.7668 0.3033 -2.528##
## Correlation matrix not shown by default, as p = 16 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it# LRT's
mm2 <- update(mm1, .~. - wbf.aud._scaled:size_pca1_scaled)
anova(mm1, mm2)## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## mm2: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +
## mm2: percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +
## mm2: stroke.amplitude_scaled:wbf.aud._scaled + stroke.amplitude_scaled:percLoad +
## mm2: stroke.amplitude_scaled:Treatment.order + stroke.amplitude_scaled:size_pca1_scaled +
## mm2: wbf.aud._scaled:percLoad + wbf.aud._scaled:Treatment.order +
## mm2: percLoad:Treatment.order + percLoad:size_pca1_scaled + Treatment.order:size_pca1_scaled
## mm1: av.resp..CO2.mL.hr. ~ (stroke.amplitude_scaled + wbf.aud._scaled +
## mm1: percLoad + Treatment.order + size_pca1_scaled)^2 + (1 | Bee.ID)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mm2 17 183.70 219.31 -74.852 149.70
## mm1 18 183.48 221.18 -73.741 147.48 2.2212 1 0.1361summary(mm2)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +
## percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +
## stroke.amplitude_scaled:wbf.aud._scaled + stroke.amplitude_scaled:percLoad +
## stroke.amplitude_scaled:Treatment.order + stroke.amplitude_scaled:size_pca1_scaled +
## wbf.aud._scaled:percLoad + wbf.aud._scaled:Treatment.order +
## percLoad:Treatment.order + percLoad:size_pca1_scaled + Treatment.order:size_pca1_scaled
## Data: bdta
##
## REML criterion at convergence: 160.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.49191 -0.41271 -0.08384 0.31098 2.15506
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 0.9109 0.9544
## Residual 0.4120 0.6419
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.3898 2.1726 0.179
## stroke.amplitude_scaled 0.9845 0.9643 1.021
## wbf.aud._scaled 0.1441 1.7214 0.084
## percLoad 6.4778 1.3979 4.634
## Treatment.order2 2.8888 2.4719 1.169
## size_pca1_scaled 3.8549 1.7749 2.172
## stroke.amplitude_scaled:wbf.aud._scaled -0.2753 0.3142 -0.876
## stroke.amplitude_scaled:percLoad -1.0519 0.6214 -1.693
## stroke.amplitude_scaled:Treatment.order2 0.1989 0.4541 0.438
## stroke.amplitude_scaled:size_pca1_scaled 0.2653 0.3414 0.777
## wbf.aud._scaled:percLoad 0.5750 1.1210 0.513
## wbf.aud._scaled:Treatment.order2 -0.3539 0.3373 -1.049
## percLoad:Treatment.order2 -1.7899 1.5426 -1.160
## percLoad:size_pca1_scaled -0.5482 1.1190 -0.490
## Treatment.order2:size_pca1_scaled -0.6698 0.2968 -2.257##
## Correlation matrix not shown by default, as p = 15 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itmm3 <- update(mm2, .~. -Treatment.order:size_pca1_scaled) # keep interaction
anova(mm3, mm2)## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## mm3: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +
## mm3: percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +
## mm3: stroke.amplitude_scaled:wbf.aud._scaled + stroke.amplitude_scaled:percLoad +
## mm3: stroke.amplitude_scaled:Treatment.order + stroke.amplitude_scaled:size_pca1_scaled +
## mm3: wbf.aud._scaled:percLoad + wbf.aud._scaled:Treatment.order +
## mm3: percLoad:Treatment.order + percLoad:size_pca1_scaled
## mm2: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +
## mm2: percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +
## mm2: stroke.amplitude_scaled:wbf.aud._scaled + stroke.amplitude_scaled:percLoad +
## mm2: stroke.amplitude_scaled:Treatment.order + stroke.amplitude_scaled:size_pca1_scaled +
## mm2: wbf.aud._scaled:percLoad + wbf.aud._scaled:Treatment.order +
## mm2: percLoad:Treatment.order + percLoad:size_pca1_scaled + Treatment.order:size_pca1_scaled
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mm3 16 188.86 222.37 -78.431 156.86
## mm2 17 183.70 219.31 -74.852 149.70 7.1571 1 0.007467 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(mm3)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +
## percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +
## stroke.amplitude_scaled:wbf.aud._scaled + stroke.amplitude_scaled:percLoad +
## stroke.amplitude_scaled:Treatment.order + stroke.amplitude_scaled:size_pca1_scaled +
## wbf.aud._scaled:percLoad + wbf.aud._scaled:Treatment.order +
## percLoad:Treatment.order + percLoad:size_pca1_scaled
## Data: bdta
##
## REML criterion at convergence: 164.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.3975 -0.5111 -0.0859 0.4055 2.1634
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 0.8748 0.9353
## Residual 0.4932 0.7023
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.4031 2.2484 0.624
## stroke.amplitude_scaled 0.4383 1.0037 0.437
## wbf.aud._scaled 0.1332 1.8467 0.072
## percLoad 5.8215 1.4465 4.025
## Treatment.order2 3.0253 2.6625 1.136
## size_pca1_scaled 3.7037 1.8855 1.964
## stroke.amplitude_scaled:wbf.aud._scaled -0.3434 0.3355 -1.024
## stroke.amplitude_scaled:percLoad -0.5762 0.6295 -0.915
## stroke.amplitude_scaled:Treatment.order2 0.1250 0.4886 0.256
## stroke.amplitude_scaled:size_pca1_scaled 0.2303 0.3602 0.639
## wbf.aud._scaled:percLoad 0.3812 1.1934 0.319
## wbf.aud._scaled:Treatment.order2 0.1745 0.2655 0.657
## percLoad:Treatment.order2 -1.8886 1.6610 -1.137
## percLoad:size_pca1_scaled -0.7152 1.1828 -0.605##
## Correlation matrix not shown by default, as p = 14 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itmm4 <- update(mm2, .~. - wbf.aud._scaled:Treatment.order)
anova(mm2, mm4)## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## mm4: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +
## mm4: percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +
## mm4: stroke.amplitude_scaled:wbf.aud._scaled + stroke.amplitude_scaled:percLoad +
## mm4: stroke.amplitude_scaled:Treatment.order + stroke.amplitude_scaled:size_pca1_scaled +
## mm4: wbf.aud._scaled:percLoad + percLoad:Treatment.order + percLoad:size_pca1_scaled +
## mm4: Treatment.order:size_pca1_scaled
## mm2: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +
## mm2: percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +
## mm2: stroke.amplitude_scaled:wbf.aud._scaled + stroke.amplitude_scaled:percLoad +
## mm2: stroke.amplitude_scaled:Treatment.order + stroke.amplitude_scaled:size_pca1_scaled +
## mm2: wbf.aud._scaled:percLoad + wbf.aud._scaled:Treatment.order +
## mm2: percLoad:Treatment.order + percLoad:size_pca1_scaled + Treatment.order:size_pca1_scaled
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mm4 16 183.54 217.05 -75.769 151.54
## mm2 17 183.70 219.31 -74.852 149.70 1.8345 1 0.1756summary(mm4)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +
## percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +
## stroke.amplitude_scaled:wbf.aud._scaled + stroke.amplitude_scaled:percLoad +
## stroke.amplitude_scaled:Treatment.order + stroke.amplitude_scaled:size_pca1_scaled +
## wbf.aud._scaled:percLoad + percLoad:Treatment.order + percLoad:size_pca1_scaled +
## Treatment.order:size_pca1_scaled
## Data: bdta
##
## REML criterion at convergence: 161.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.54564 -0.45657 -0.08531 0.37288 2.24119
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 0.8871 0.9418
## Residual 0.4227 0.6502
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.09727 2.17454 0.045
## stroke.amplitude_scaled 0.67101 0.92807 0.723
## wbf.aud._scaled 0.01824 1.73163 0.011
## percLoad 6.66930 1.39819 4.770
## Treatment.order2 3.77686 2.36346 1.598
## size_pca1_scaled 3.69859 1.77947 2.078
## stroke.amplitude_scaled:wbf.aud._scaled -0.29890 0.31607 -0.946
## stroke.amplitude_scaled:percLoad -0.88041 0.60607 -1.453
## stroke.amplitude_scaled:Treatment.order2 0.29973 0.44709 0.670
## stroke.amplitude_scaled:size_pca1_scaled 0.19508 0.33683 0.579
## wbf.aud._scaled:percLoad 0.52625 1.12791 0.467
## percLoad:Treatment.order2 -2.35292 1.47151 -1.599
## percLoad:size_pca1_scaled -0.53345 1.12492 -0.474
## Treatment.order2:size_pca1_scaled -0.45825 0.21875 -2.095##
## Correlation matrix not shown by default, as p = 14 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itmm5 <- update(mm4, .~. - stroke.amplitude_scaled:wbf.aud._scaled)
anova(mm4, mm5)## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## mm5: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +
## mm5: percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +
## mm5: stroke.amplitude_scaled:percLoad + stroke.amplitude_scaled:Treatment.order +
## mm5: stroke.amplitude_scaled:size_pca1_scaled + wbf.aud._scaled:percLoad +
## mm5: percLoad:Treatment.order + percLoad:size_pca1_scaled + Treatment.order:size_pca1_scaled
## mm4: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +
## mm4: percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +
## mm4: stroke.amplitude_scaled:wbf.aud._scaled + stroke.amplitude_scaled:percLoad +
## mm4: stroke.amplitude_scaled:Treatment.order + stroke.amplitude_scaled:size_pca1_scaled +
## mm4: wbf.aud._scaled:percLoad + percLoad:Treatment.order + percLoad:size_pca1_scaled +
## mm4: Treatment.order:size_pca1_scaled
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mm5 15 182.56 213.98 -76.281 152.56
## mm4 16 183.54 217.05 -75.769 151.54 1.0242 1 0.3115summary(mm5)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +
## percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +
## stroke.amplitude_scaled:percLoad + stroke.amplitude_scaled:Treatment.order +
## stroke.amplitude_scaled:size_pca1_scaled + wbf.aud._scaled:percLoad +
## percLoad:Treatment.order + percLoad:size_pca1_scaled + Treatment.order:size_pca1_scaled
## Data: bdta
##
## REML criterion at convergence: 161.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6390 -0.4133 -0.1268 0.3654 2.3371
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 0.9212 0.9598
## Residual 0.4080 0.6388
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.1984 2.1527 0.092
## stroke.amplitude_scaled 0.9261 0.8850 1.047
## wbf.aud._scaled 1.4238 0.8461 1.683
## percLoad 6.6563 1.3848 4.807
## Treatment.order2 4.0182 2.3080 1.741
## size_pca1_scaled 4.4330 1.5838 2.799
## stroke.amplitude_scaled:percLoad -1.0551 0.5757 -1.833
## stroke.amplitude_scaled:Treatment.order2 0.3835 0.4335 0.885
## stroke.amplitude_scaled:size_pca1_scaled 0.3313 0.3041 1.089
## wbf.aud._scaled:percLoad -0.3826 0.5739 -0.667
## percLoad:Treatment.order2 -2.5079 1.4368 -1.745
## percLoad:size_pca1_scaled -1.0040 1.0016 -1.002
## Treatment.order2:size_pca1_scaled -0.4638 0.2153 -2.154##
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itmm6 <- update(mm5, .~. - stroke.amplitude_scaled:Treatment.order)
anova(mm5, mm6)## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## mm6: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +
## mm6: percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +
## mm6: stroke.amplitude_scaled:percLoad + stroke.amplitude_scaled:size_pca1_scaled +
## mm6: wbf.aud._scaled:percLoad + percLoad:Treatment.order + percLoad:size_pca1_scaled +
## mm6: Treatment.order:size_pca1_scaled
## mm5: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +
## mm5: percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +
## mm5: stroke.amplitude_scaled:percLoad + stroke.amplitude_scaled:Treatment.order +
## mm5: stroke.amplitude_scaled:size_pca1_scaled + wbf.aud._scaled:percLoad +
## mm5: percLoad:Treatment.order + percLoad:size_pca1_scaled + Treatment.order:size_pca1_scaled
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mm6 14 181.81 211.13 -76.903 153.81
## mm5 15 182.56 213.98 -76.281 152.56 1.2439 1 0.2647summary(mm6)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +
## percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +
## stroke.amplitude_scaled:percLoad + stroke.amplitude_scaled:size_pca1_scaled +
## wbf.aud._scaled:percLoad + percLoad:Treatment.order + percLoad:size_pca1_scaled +
## Treatment.order:size_pca1_scaled
## Data: bdta
##
## REML criterion at convergence: 162.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.67475 -0.46717 -0.06738 0.39201 2.38883
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 0.8984 0.9479
## Residual 0.4128 0.6425
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.9357 1.9719 0.474
## stroke.amplitude_scaled 1.2688 0.7898 1.607
## wbf.aud._scaled 1.5030 0.8458 1.777
## percLoad 6.2134 1.2846 4.837
## Treatment.order2 2.7598 1.7714 1.558
## size_pca1_scaled 4.3302 1.5810 2.739
## stroke.amplitude_scaled:percLoad -1.1420 0.5668 -2.015
## stroke.amplitude_scaled:size_pca1_scaled 0.2958 0.3017 0.981
## wbf.aud._scaled:percLoad -0.4589 0.5695 -0.806
## percLoad:Treatment.order2 -1.7307 1.1093 -1.560
## percLoad:size_pca1_scaled -0.9794 1.0021 -0.977
## Treatment.order2:size_pca1_scaled -0.4111 0.2072 -1.984
##
## Correlation of Fixed Effects:
## (Intr) strk._ wbf.._ percLd Trtm.2 sz_p1_ st._:L s._:_1 w.._:L
## strk.mpltd_ -0.223
## wbf.d._scld -0.060 -0.037
## percLoad -0.991 0.283 0.065
## Trtmnt.rdr2 -0.545 0.057 0.375 0.530
## sz_pc1_scld -0.191 0.198 0.441 0.174 0.153
## strk.mpl_:L 0.485 -0.939 0.062 -0.542 -0.082 -0.207
## strk.m_:_1_ 0.005 0.214 0.080 -0.022 -0.105 0.841 -0.174
## wbf.d._sc:L 0.161 -0.087 -0.967 -0.178 -0.396 -0.427 0.075 -0.072
## prcLd:Trt.2 0.562 -0.070 -0.370 -0.554 -0.994 -0.132 0.101 0.125 0.405
## prcLd:sz_1_ 0.170 -0.236 -0.461 -0.158 -0.138 -0.985 0.231 -0.852 0.469
## Trtmn.2:_1_ 0.378 -0.161 0.198 -0.381 -0.218 0.094 0.258 0.240 -0.190
## pL:T.2 pL:_1_
## strk.mpltd_
## wbf.d._scld
## percLoad
## Trtmnt.rdr2
## sz_pc1_scld
## strk.mpl_:L
## strk.m_:_1_
## wbf.d._sc:L
## prcLd:Trt.2
## prcLd:sz_1_ 0.122
## Trtmn.2:_1_ 0.219 -0.171mm7 <- update(mm6, .~. - percLoad:Treatment.order)
anova(mm6, mm7)## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## mm7: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +
## mm7: percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +
## mm7: stroke.amplitude_scaled:percLoad + stroke.amplitude_scaled:size_pca1_scaled +
## mm7: wbf.aud._scaled:percLoad + percLoad:size_pca1_scaled + Treatment.order:size_pca1_scaled
## mm6: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +
## mm6: percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +
## mm6: stroke.amplitude_scaled:percLoad + stroke.amplitude_scaled:size_pca1_scaled +
## mm6: wbf.aud._scaled:percLoad + percLoad:Treatment.order + percLoad:size_pca1_scaled +
## mm6: Treatment.order:size_pca1_scaled
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mm7 13 181.80 209.02 -77.899 155.80
## mm6 14 181.81 211.13 -76.903 153.81 1.9911 1 0.1582summary(mm7)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +
## percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +
## stroke.amplitude_scaled:percLoad + stroke.amplitude_scaled:size_pca1_scaled +
## wbf.aud._scaled:percLoad + percLoad:size_pca1_scaled + Treatment.order:size_pca1_scaled
## Data: bdta
##
## REML criterion at convergence: 167
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.71405 -0.37038 -0.02943 0.33275 2.39389
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.0729 1.0358
## Residual 0.3714 0.6094
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.634514 1.608837 1.638
## stroke.amplitude_scaled 1.315040 0.769176 1.710
## wbf.aud._scaled 1.049892 0.754406 1.392
## percLoad 5.141776 1.055823 4.870
## Treatment.order2 -0.002183 0.188595 -0.012
## size_pca1_scaled 4.065746 1.542299 2.636
## stroke.amplitude_scaled:percLoad -1.137541 0.551286 -2.063
## stroke.amplitude_scaled:size_pca1_scaled 0.361927 0.296016 1.223
## wbf.aud._scaled:percLoad -0.150894 0.503708 -0.300
## percLoad:size_pca1_scaled -0.850048 0.981109 -0.866
## Treatment.order2:size_pca1_scaled -0.336078 0.193590 -1.736
##
## Correlation of Fixed Effects:
## (Intr) strk._ wbf.._ percLd Trtm.2 sz_p1_ st._:L s._:_1 w.._:L
## strk.mpltd_ -0.226
## wbf.d._scld 0.179 -0.049
## percLoad -0.987 0.297 -0.168
## Trtmnt.rdr2 0.165 -0.132 0.069 -0.237
## sz_pc1_scld -0.126 0.212 0.432 0.107 0.194
## strk.mpl_:L 0.528 -0.938 0.089 -0.593 0.188 -0.207
## strk.m_:_1_ -0.053 0.242 0.152 0.032 0.169 0.877 -0.192
## wbf.d._sc:L -0.067 -0.094 -0.959 0.039 0.072 -0.421 0.068 -0.154
## prcLd:sz_1_ 0.109 -0.253 -0.458 -0.097 -0.151 -0.984 0.235 -0.887 0.472
## Trtmn.2:_1_ 0.324 -0.134 0.311 -0.327 -0.006 0.145 0.235 0.239 -0.314
## pL:_1_
## strk.mpltd_
## wbf.d._scld
## percLoad
## Trtmnt.rdr2
## sz_pc1_scld
## strk.mpl_:L
## strk.m_:_1_
## wbf.d._sc:L
## prcLd:sz_1_
## Trtmn.2:_1_ -0.221mm8 <- update(mm7, .~. - wbf.aud._scaled:percLoad)
anova(mm8, mm7)## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## mm8: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +
## mm8: percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +
## mm8: stroke.amplitude_scaled:percLoad + stroke.amplitude_scaled:size_pca1_scaled +
## mm8: percLoad:size_pca1_scaled + Treatment.order:size_pca1_scaled
## mm7: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +
## mm7: percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +
## mm7: stroke.amplitude_scaled:percLoad + stroke.amplitude_scaled:size_pca1_scaled +
## mm7: wbf.aud._scaled:percLoad + percLoad:size_pca1_scaled + Treatment.order:size_pca1_scaled
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mm8 12 179.99 205.12 -77.996 155.99
## mm7 13 181.80 209.02 -77.899 155.80 0.1943 1 0.6594summary(mm8)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +
## percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +
## stroke.amplitude_scaled:percLoad + stroke.amplitude_scaled:size_pca1_scaled +
## percLoad:size_pca1_scaled + Treatment.order:size_pca1_scaled
## Data: bdta
##
## REML criterion at convergence: 167.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.73514 -0.34367 -0.05258 0.32887 2.40909
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.0586 1.0289
## Residual 0.3631 0.6026
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.600949 1.589185 1.637
## stroke.amplitude_scaled 1.297131 0.757782 1.712
## wbf.aud._scaled 0.831911 0.211132 3.940
## percLoad 5.155421 1.044492 4.936
## Treatment.order2 0.001467 0.186058 0.008
## size_pca1_scaled 3.871007 1.384576 2.796
## stroke.amplitude_scaled:percLoad -1.128710 0.544304 -2.074
## stroke.amplitude_scaled:size_pca1_scaled 0.348367 0.289558 1.203
## percLoad:size_pca1_scaled -0.711749 0.856171 -0.831
## Treatment.order2:size_pca1_scaled -0.354355 0.181766 -1.950
##
## Correlation of Fixed Effects:
## (Intr) strk._ wbf.._ percLd Trtm.2 sz_p1_ st._:L s._:_1 pL:_1_
## strk.mpltd_ -0.234
## wbf.d._scld 0.407 -0.497
## percLoad -0.987 0.302 -0.462
## Trtmnt.rdr2 0.172 -0.126 0.492 -0.241
## sz_pc1_scld -0.169 0.191 0.107 0.136 0.248
## strk.mpl_:L 0.535 -0.938 0.549 -0.598 0.184 -0.196
## strk.m_:_1_ -0.063 0.231 0.015 0.037 0.183 0.906 -0.184
## prcLd:sz_1_ 0.160 -0.238 -0.020 -0.131 -0.211 -0.982 0.231 -0.935
## Trtmn.2:_1_ 0.320 -0.173 0.034 -0.332 0.018 0.015 0.270 0.204 -0.087mm9 <- update(mm8, .~. - percLoad:size_pca1_scaled)
anova(mm9, mm8)## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## mm9: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +
## mm9: percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +
## mm9: stroke.amplitude_scaled:percLoad + stroke.amplitude_scaled:size_pca1_scaled +
## mm9: Treatment.order:size_pca1_scaled
## mm8: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +
## mm8: percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +
## mm8: stroke.amplitude_scaled:percLoad + stroke.amplitude_scaled:size_pca1_scaled +
## mm8: percLoad:size_pca1_scaled + Treatment.order:size_pca1_scaled
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mm9 11 178.87 201.91 -78.436 156.87
## mm8 12 179.99 205.12 -77.996 155.99 0.8806 1 0.348summary(mm9)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +
## percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +
## stroke.amplitude_scaled:percLoad + stroke.amplitude_scaled:size_pca1_scaled +
## Treatment.order:size_pca1_scaled
## Data: bdta
##
## REML criterion at convergence: 169.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.87306 -0.37383 -0.07975 0.32076 2.55486
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.0512 1.025
## Residual 0.3612 0.601
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.81219 1.56426 1.798
## stroke.amplitude_scaled 1.14638 0.73391 1.562
## wbf.aud._scaled 0.82870 0.21048 3.937
## percLoad 5.04157 1.03253 4.883
## Treatment.order2 -0.03107 0.18138 -0.171
## size_pca1_scaled 2.74118 0.26284 10.429
## stroke.amplitude_scaled:percLoad -1.02381 0.52811 -1.939
## stroke.amplitude_scaled:size_pca1_scaled 0.12336 0.10254 1.203
## Treatment.order2:size_pca1_scaled -0.36752 0.18058 -2.035
##
## Correlation of Fixed Effects:
## (Intr) strk._ wbf.._ percLd Trtm.2 sz_p1_ st._:L s._:_1
## strk.mpltd_ -0.204
## wbf.d._scld 0.415 -0.516
## percLoad -0.987 0.282 -0.468
## Trtmnt.rdr2 0.213 -0.186 0.499 -0.277
## sz_pc1_scld -0.068 -0.234 0.463 0.039 0.222
## strk.mpl_:L 0.518 -0.934 0.569 -0.589 0.245 0.163
## strk.m_:_1_ 0.245 0.025 -0.008 -0.241 -0.042 -0.170 0.091
## Trtmn.2:_1_ 0.340 -0.200 0.033 -0.348 0.000 -0.374 0.300 0.345mm10 <- update(mm9, .~. - stroke.amplitude_scaled:percLoad)
anova(mm10, mm9) # keep interaction## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## mm10: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +
## mm10: percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +
## mm10: stroke.amplitude_scaled:size_pca1_scaled + Treatment.order:size_pca1_scaled
## mm9: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +
## mm9: percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +
## mm9: stroke.amplitude_scaled:percLoad + stroke.amplitude_scaled:size_pca1_scaled +
## mm9: Treatment.order:size_pca1_scaled
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mm10 10 181.17 202.12 -80.587 161.17
## mm9 11 178.87 201.91 -78.436 156.87 4.3005 1 0.0381 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(mm10)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +
## percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +
## stroke.amplitude_scaled:size_pca1_scaled + Treatment.order:size_pca1_scaled
## Data: bdta
##
## REML criterion at convergence: 173.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.92324 -0.37734 -0.04338 0.31256 2.70415
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 0.9448 0.9720
## Residual 0.4257 0.6524
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 4.27791 1.40239 3.050
## stroke.amplitude_scaled -0.20299 0.27398 -0.741
## wbf.aud._scaled 1.06817 0.18262 5.849
## percLoad 3.92936 0.87680 4.481
## Treatment.order2 0.05687 0.18987 0.300
## size_pca1_scaled 2.83452 0.26014 10.896
## stroke.amplitude_scaled:size_pca1_scaled 0.14257 0.10966 1.300
## Treatment.order2:size_pca1_scaled -0.26620 0.18609 -1.431
##
## Correlation of Fixed Effects:
## (Intr) strk._ wbf.._ percLd Trtm.2 sz_p1_ s._:_1
## strk.mpltd_ 0.918
## wbf.d._scld 0.183 0.066
## percLoad -0.988 -0.927 -0.213
## Trtmnt.rdr2 0.102 0.129 0.442 -0.170
## sz_pc1_scld -0.182 -0.233 0.476 0.169 0.195
## strk.m_:_1_ 0.214 0.291 -0.061 -0.214 -0.064 -0.189
## Trtmn.2:_1_ 0.215 0.222 -0.167 -0.211 -0.076 -0.474 0.326mm11 <- update(mm9, .~. - stroke.amplitude_scaled:size_pca1_scaled)
anova(mm9, mm11) # drop interaction## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## mm11: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +
## mm11: percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +
## mm11: stroke.amplitude_scaled:percLoad + Treatment.order:size_pca1_scaled
## mm9: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +
## mm9: percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +
## mm9: stroke.amplitude_scaled:percLoad + stroke.amplitude_scaled:size_pca1_scaled +
## mm9: Treatment.order:size_pca1_scaled
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mm11 10 178.60 199.54 -79.300 158.60
## mm9 11 178.87 201.91 -78.436 156.87 1.7279 1 0.1887summary(mm11)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +
## percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +
## stroke.amplitude_scaled:percLoad + Treatment.order:size_pca1_scaled
## Data: bdta
##
## REML criterion at convergence: 168.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.75223 -0.41953 -0.02529 0.44228 2.43853
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.0672 1.0331
## Residual 0.3623 0.6019
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.34652 1.52140 1.542
## stroke.amplitude_scaled 1.12885 0.73555 1.535
## wbf.aud._scaled 0.82929 0.21120 3.927
## percLoad 5.34378 1.00541 5.315
## Treatment.order2 -0.02245 0.18159 -0.124
## size_pca1_scaled 2.79423 0.26035 10.732
## stroke.amplitude_scaled:percLoad -1.08484 0.52736 -2.057
## Treatment.order2:size_pca1_scaled -0.44257 0.16981 -2.606
##
## Correlation of Fixed Effects:
## (Intr) strk._ wbf.._ percLd Trtm.2 sz_p1_ st._:L
## strk.mpltd_ -0.217
## wbf.d._scld 0.431 -0.517
## percLoad -0.986 0.297 -0.485
## Trtmnt.rdr2 0.231 -0.186 0.500 -0.296
## sz_pc1_scld -0.027 -0.233 0.468 -0.002 0.218
## strk.mpl_:L 0.514 -0.941 0.573 -0.587 0.251 0.182
## Trtmn.2:_1_ 0.281 -0.222 0.038 -0.291 0.016 -0.339 0.287## main effects
mm14 <- update(mm11, .~. - wbf.aud._scaled)
anova(mm11, mm14) # wbf is significant## refitting model(s) with ML (instead of REML)## Data: bdta
## Models:
## mm14: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + percLoad + Treatment.order +
## mm14: size_pca1_scaled + (1 | Bee.ID) + stroke.amplitude_scaled:percLoad +
## mm14: Treatment.order:size_pca1_scaled
## mm11: av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +
## mm11: percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +
## mm11: stroke.amplitude_scaled:percLoad + Treatment.order:size_pca1_scaled
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mm14 9 191.71 210.56 -86.855 173.71
## mm11 10 178.60 199.54 -79.300 158.60 15.11 1 0.0001014 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(mm11)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +
## percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +
## stroke.amplitude_scaled:percLoad + Treatment.order:size_pca1_scaled
## Data: bdta
##
## REML criterion at convergence: 168.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.75223 -0.41953 -0.02529 0.44228 2.43853
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.0672 1.0331
## Residual 0.3623 0.6019
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.34652 1.52140 1.542
## stroke.amplitude_scaled 1.12885 0.73555 1.535
## wbf.aud._scaled 0.82929 0.21120 3.927
## percLoad 5.34378 1.00541 5.315
## Treatment.order2 -0.02245 0.18159 -0.124
## size_pca1_scaled 2.79423 0.26035 10.732
## stroke.amplitude_scaled:percLoad -1.08484 0.52736 -2.057
## Treatment.order2:size_pca1_scaled -0.44257 0.16981 -2.606
##
## Correlation of Fixed Effects:
## (Intr) strk._ wbf.._ percLd Trtm.2 sz_p1_ st._:L
## strk.mpltd_ -0.217
## wbf.d._scld 0.431 -0.517
## percLoad -0.986 0.297 -0.485
## Trtmnt.rdr2 0.231 -0.186 0.500 -0.296
## sz_pc1_scld -0.027 -0.233 0.468 -0.002 0.218
## strk.mpl_:L 0.514 -0.941 0.573 -0.587 0.251 0.182
## Trtmn.2:_1_ 0.281 -0.222 0.038 -0.291 0.016 -0.339 0.287summary(mm11) # final model for paper## Linear mixed model fit by REML ['lmerMod']
## Formula:
## av.resp..CO2.mL.hr. ~ stroke.amplitude_scaled + wbf.aud._scaled +
## percLoad + Treatment.order + size_pca1_scaled + (1 | Bee.ID) +
## stroke.amplitude_scaled:percLoad + Treatment.order:size_pca1_scaled
## Data: bdta
##
## REML criterion at convergence: 168.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.75223 -0.41953 -0.02529 0.44228 2.43853
##
## Random effects:
## Groups Name Variance Std.Dev.
## Bee.ID (Intercept) 1.0672 1.0331
## Residual 0.3623 0.6019
## Number of obs: 60, groups: Bee.ID, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.34652 1.52140 1.542
## stroke.amplitude_scaled 1.12885 0.73555 1.535
## wbf.aud._scaled 0.82929 0.21120 3.927
## percLoad 5.34378 1.00541 5.315
## Treatment.order2 -0.02245 0.18159 -0.124
## size_pca1_scaled 2.79423 0.26035 10.732
## stroke.amplitude_scaled:percLoad -1.08484 0.52736 -2.057
## Treatment.order2:size_pca1_scaled -0.44257 0.16981 -2.606
##
## Correlation of Fixed Effects:
## (Intr) strk._ wbf.._ percLd Trtm.2 sz_p1_ st._:L
## strk.mpltd_ -0.217
## wbf.d._scld 0.431 -0.517
## percLoad -0.986 0.297 -0.485
## Trtmnt.rdr2 0.231 -0.186 0.500 -0.296
## sz_pc1_scld -0.027 -0.233 0.468 -0.002 0.218
## strk.mpl_:L 0.514 -0.941 0.573 -0.587 0.251 0.182
## Trtmn.2:_1_ 0.281 -0.222 0.038 -0.291 0.016 -0.339 0.287# write output
summary(mm11)$coefficients ## Estimate Std. Error t value
## (Intercept) 2.34652139 1.5213993 1.5423442
## stroke.amplitude_scaled 1.12884653 0.7355504 1.5346963
## wbf.aud._scaled 0.82929053 0.2111955 3.9266482
## percLoad 5.34378104 1.0054092 5.3150309
## Treatment.order2 -0.02245326 0.1815899 -0.1236482
## size_pca1_scaled 2.79422560 0.2603545 10.7323894
## stroke.amplitude_scaled:percLoad -1.08484103 0.5273643 -2.0570998
## Treatment.order2:size_pca1_scaled -0.44257485 0.1698069 -2.6063418write.csv( summary(mm11)$coefficients, file = "resp_oth_Coefs_percLoad.csv" )
# diagnostics
# qq plot
qqnorm(resid(mm11), main = "")
qqline(resid(mm11)) # good, but two outliers
# residual plot
plot(fitted(mm11), resid(mm11), xlab = "fitted", ylab = "residuals")
abline(0,0)
# QQPlot for group-level effects
qqnorm(ranef(mm11)$Bee.ID[[1]], main="Normal Q-Q plot for random effects")
qqline(ranef(mm11)$Bee.ID[[1]]) # looks good
This model is fundamentally different from the previous one. It’s saying, when we account for differences in wingbeat freq, stroke amplitude, and bee size, we sill see an effect of percent load on metabolic rate.
OR, we could say that holding all other variables constant, an increase in wingbeat frequency is associated with an increase in metabolic rate.